//-----------------------------------------------------------------------------
//
// Date : 24/05/02
// Author : Luis M. de la Cruz.
// Description: Conveccion-Difusion en 2D completa
//

#include <cmath>
#include <algorithm>

#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS_CoDi.hpp"
#include "Equations/Momentum_XCoDi.hpp"
#include "Schemes/CDS_XCoDi.hpp"
#include "Equations/Momentum_YCoDi.hpp"
#include "Schemes/CDS_YCoDi.hpp"
#include "Equations/PressureCorrection.hpp"
#include "Schemes/Simplec.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
#include "Utils/diff_operators.hpp"
using namespace Tuna;

#define PI 3.14159265358979323846

// Scalar field in two dimensions
typedef TunaArray<double,2>::huge ScalarField2D;

//
// DT = 0.0001, a complete cicle is 800
// 
//int intervalo_1 = 50, intervalo_2 = 300;
//
// DT = 0.00001, a complete cicle is 8000, imon = 64, jmon = 192
// 257 x 257
//
// int intervalo_1 = 500, intervalo_2 = 3000;
//
// DT = 0.00001, , a complete cicle is 80000
//
//int intervalo_1 = 5000, intervalo_2 = 30000;

int intervalo_1, intervalo_2;
double freq;
int ciclo;


inline void boundaryConditions(ScalarField2D &T, int iteracion) {

    int hot_up, hot_flat, hot_down, col_up, col_flat, col_down;
    int i, opcion;
    int it1 = intervalo_1;
    int it2 = it1 + intervalo_2;
    int it3 = it2 + intervalo_1;
    int it4 = it3 + intervalo_1;
    int it5 = it4 + intervalo_2;
    int it6 = it5 + intervalo_1;
    double thot, tcol;

    hot_up = hot_flat = hot_down = col_up = col_flat = col_down = 0;
    opcion = iteracion % ciclo;

    if (0 < opcion and opcion <= it1) {
	hot_up = 1;
	i = opcion;
    } 
    else if (it1 < opcion and opcion <= it2) {
	hot_flat = 1;
	i = opcion - it1;
    }
    else if (it2 < opcion and opcion <= it3) {
	hot_down = 1;
	i = opcion - it2;
    }
    else if (it3 < opcion and opcion <= it4) {
	col_up = 1;
	i = opcion - it3;
    }
    else if (it4 < opcion and opcion <= it5) {
	col_flat = 1;
	i = opcion - it4;
    }
    else if (it5 < opcion and opcion <= it6 or opcion == 0) {
	col_down = 1;
	i = opcion - it5;
    }

    
    thot =  0.5 * ( hot_up * sin(freq*i) * sin(freq*i) +
		    hot_flat +
		    hot_down * sin(freq*(intervalo_1+i)) * 
		               sin(freq*(intervalo_1+i)) );

    tcol = -0.5 * ( col_up * sin(freq*i) * sin(freq * i) +
		    col_flat +
		    col_down * sin(freq*(intervalo_1+i)) * 
		               sin(freq*(intervalo_1+i)) );    


    /***
    int num_vols_x = T.ubound(firstDim) + 1;
    Range I1(0,num_vols_x / 2 - 1);
    Range I2(num_vols_x / 2, num_vols_x-1);
    T(I2, T.lbound(secondDim)) = thot;  
    T(I1, T.ubound(secondDim)) = tcol;
    ***/

    Range all = Range::all();
    T(T.lbound(firstDim), all) = tcol;  
    T(T.ubound(firstDim), all) = thot; 
}


template<class T_energy, class T_xmom, class T_ymom, class T_press>
double SIMPLEC(T_energy &, T_xmom &, T_ymom &, T_press &, 
	       int, double);

double errorT, errorX, errorY, errorP;

int stop;

int main() 
{

    double lambda, Amplitud;
    double top_wall, bot_wall;
    double longitud_x, longitud_y;
    int num_nodos_x, num_nodos_y;
    double tolerancia, dt;    
    int max_iter, iteraciones_max, frecuencia;
    double alpha_p, alpha_u, alpha_v;
    double Prandtl, Rayleigh;

    int imon, jmon;

// Reading data from the input file
    std::ifstream input_cin ("input7");
    input_cin >> longitud_x
	      >> longitud_y
	      >> num_nodos_x 
	      >> num_nodos_y
	      >> top_wall
	      >> bot_wall 
	      >> dt
	      >> iteraciones_max
	      >> tolerancia
	      >> max_iter
	      >> Prandtl
	      >> Rayleigh
	      >> frecuencia
	      >> imon
	      >> jmon
	      >> intervalo_1 
	      >> intervalo_2;

    freq = PI / (2 * intervalo_1); 
    ciclo = 4 * intervalo_1 + 2 * intervalo_2;

//
//  Linear System 
//
    SparseMatrix< Diagonal<double, 2> > A(num_nodos_x, num_nodos_y); 

    ScalarField2D             b(num_nodos_x, num_nodos_y);
//
//  Mesh definition
//
    StructuredMesh<Uniform<double, 2> > malla(longitud_x, num_nodos_x,
					      longitud_y, num_nodos_y);
    malla.print();
    double dx = malla.getDelta(X);    
    double dy = malla.getDelta(Y);
    int num_vols_x = malla.getExtentVolumes(X);
    int num_vols_y = malla.getExtentVolumes(Y);  
//
//  Scalar fields definitions (in volumes centers)
//
    ScalarField2D T(malla.getExtentVolumes());      // Temperature
    ScalarField2D p(malla.getExtentVolumes());      // pressure
    ScalarField2D pp(malla.getExtentVolumes());     // pressure correction
    ScalarField2D u(malla.getExtentVolumes());      // u-velocity
    ScalarField2D v(malla.getExtentVolumes());      // v-velocity
    ScalarField2D vort(malla.getExtentVolumes());   // vorticity
    ScalarField2D us(num_nodos_x, num_vols_y); // u-velocity staggered
    ScalarField2D vs(num_vols_x, num_nodos_y); // v-velocity staggered
//
//  Scalar fields definitions for printing (in nodes)
//
    ScalarField2D Tn(malla.getExtentNodes());      // Temperature
    ScalarField2D pn(malla.getExtentNodes());      // pressure
    ScalarField2D un(malla.getExtentNodes());      // u-velocity
    ScalarField2D vn(malla.getExtentNodes());      // v-velocity
    ScalarField2D vortn(malla.getExtentNodes());   // vorticity
//
// Energy equation definition
//
    ScalarEquation<CDS_CoDi<double,2> > energia(T, A, b, malla.getDeltas());
    energia.setDeltaTime(dt);
    energia.setGamma(1.0);
    energia.setDirichlet(LEFT_WALL);
    energia.setDirichlet(RIGHT_WALL);
    energia.setNeumann(TOP_WALL);
    energia.setNeumann(BOTTOM_WALL);
    energia.setUvelocity(us);
    energia.setVvelocity(vs);
    energia.print();
//
// Momentum-x equation definition
//
    Momentum_XCoDi<CDS_XCoDi<double, 2> > momento_x(us, A, b, malla.getDeltas());
    momento_x.setGamma(Prandtl);
    momento_x.setDeltaTime(dt);
    momento_x.setDirichlet(LEFT_WALL, 0);
    momento_x.setDirichlet(RIGHT_WALL, 0);
    momento_x.setDirichlet(TOP_WALL, 0);
    momento_x.setDirichlet(BOTTOM_WALL, 0);
    momento_x.setUvelocity(us);
    momento_x.setVvelocity(vs);
    momento_x.setPressure(p);
    momento_x.print();
//    
// Momentum-x equation definition
//
    Momentum_YCoDi<CDS_YCoDi<double, 2> > momento_y(vs, A, b, malla.getDeltas());
    momento_y.setGamma(Prandtl);
    momento_y.setDeltaTime(dt);
    momento_y.setDirichlet(LEFT_WALL, 0);
    momento_y.setDirichlet(RIGHT_WALL, 0);
    momento_y.setDirichlet(TOP_WALL, 0);
    momento_y.setDirichlet(BOTTOM_WALL, 0);
    momento_y.setUvelocity(us);
    momento_y.setVvelocity(vs);
    momento_y.setPressure(p);
    momento_y.setTemperature(T);
    momento_y.setRayleigh(Rayleigh);
    momento_y.print();
//
// Pressure equation definition
//
    PressureCorrection<Simplec<double, 2 > > presion(pp, A, b, malla.getDeltas());
    presion.setDeltaTime(dt);
    presion.applyBounds(1, num_nodos_x-1, 1, num_nodos_y-1);
    presion.setUvelocity(us);
    presion.setVvelocity(vs);
    presion.set_du(momento_x.get_du());
    presion.set_dv(momento_y.get_dv());
    presion.setPressure(p);
    presion.print();

    ofstream sorsumfile("./Data7/sorsum"), 
	nusselt_hot_file("./Data7/nuh"), nusselt_cold_file("./Data7/nuc"),
	umon_file("./Data7/umon"), vmon_file("./Data7/vmon"), 
	Tmon_file("./Data7/tmon"), pmon_file("./Data7/pmon"),
	errorTfile("./Data7/errorT"), errorPfile("./Data7/errorP"),
	errorXfile("./Data7/errorX"), errorYfile("./Data7/errorY");

    double sorsum, resorm, nuh, nuc;
    
    int iteracion = 0;

    NumUtils::interpolateToNodes(Tn, T);
    NumUtils::interpolateToNodes(pn, p);
    NumUtils::interpolateToNodesUs(un, us);
    NumUtils::interpolateToNodesVs(vn, vs);
    InOut::writeToFile_DX(Tn, iteracion , "./Data7/temp.", dx, dy);
    InOut::writeToFile_DX(pn, iteracion , "./Data7/pres.", dx, dy);
    InOut::writeToFile_DX(un, vn, iteracion , "./Data7/velc.", dx, dy);
    

    for(iteracion = 1; iteracion <= iteraciones_max; ++iteracion) {

        boundaryConditions(T, iteracion);
	energia.updatePhi(T);
      
        sorsum = SIMPLEC(energia, momento_x, momento_y, presion,
			 max_iter, tolerancia);
	

	if ( !(iteracion % frecuencia) ) {
	    NumUtils::interpolateToNodes(Tn, T);
	    NumUtils::interpolateToNodes(pn, p);
	    NumUtils::interpolateToNodesUs(un, us);
	    NumUtils::interpolateToNodesVs(vn, vs);

	    InOut::writeToFile_DX(Tn, iteracion , "./Data7/temp.", dx, dy);
	    InOut::writeToFile_DX(pn, iteracion, "./Data7/pres.", dx, dy);
	    InOut::writeToFile_DX(un, vn, iteracion, "./Data7/velc.", dx, dy);

	    DiffOperators::vorticity(vortn, un, vn, dx, dy);
	    InOut::writeToFile_DX(vortn, iteracion, "./Data7/vort.", dx, dy);
	}

	if(!(iteracion % 1)) { 
	    sorsumfile << iteracion << "\t" << sorsum << endl;
	    
	    nuh = NumUtils::calcNusseltHoriz(T, dx, dy, 0);
	    nuc = NumUtils::calcNusseltHoriz(T, dx, dy, num_vols_y - 2);
	    
	    nusselt_hot_file << iteracion << "\t" << nuh << endl;
	    nusselt_cold_file << iteracion << "\t" << nuc << endl;
	    errorTfile << iteracion << "\t" << errorT << endl;
	    errorPfile << iteracion << "\t" << errorP << endl;
	    errorXfile << iteracion << "\t" << errorX << endl;
	    errorYfile << iteracion << "\t" << errorY << endl;
	    
	    umon_file << iteracion << "\t" << un(imon,jmon) << endl;
	    vmon_file << iteracion << "\t" << vn(imon,jmon) << endl;
	    Tmon_file << iteracion << "\t" << Tn(imon,jmon) << endl;
	    pmon_file << iteracion << "\t" << pn(imon,jmon) << endl;
	    
	    cout << "=> Iter = " << iteracion 
		 << " s = " << sorsum 
		 << " nuc = " << nuc 
		 << " nuh = " << nuh << endl
		 << " e_T = " << errorT
		 << " e_P = " << errorP
		 << " e_X = " << errorX
		 << " e_Y = " << errorY << endl;
	}
    }

    std::cout << "\n\n NOTE: The data were written in the [./Data7] directory in DX format";
    return 0;
}

template<class T_energy, class T_xmom, class T_ymom, class T_press>
double SIMPLEC(T_energy &energy, T_xmom &xmomentum, 
	       T_ymom &ymomentum, T_press &pressure, 
	       int max_iter, double tolerancia)
{
    double sorsum = 10.0, tolerance = 1e-2;
    int counter = 0;

    while ( (sorsum > tolerance) && (counter < 20) ) {

	energy.calcCoefficients();
	Solver::TDMA2DY(energy, tolerancia, max_iter);
	errorT = energy.calcErrorL1();
	energy.update();

	xmomentum.calcCoefficients();
	Solver::TDMA2DY(xmomentum, tolerancia, max_iter);
	errorX = xmomentum.calcErrorL1();
	xmomentum.update();

	ymomentum.calcCoefficients();
	Solver::TDMA2DY(ymomentum, tolerancia, max_iter);
	errorY = ymomentum.calcErrorL1();
	ymomentum.update();

	pressure.calcCoefficients();
	Solver::TDMA2DY(pressure, tolerancia, max_iter);
	pressure.correction();	
	sorsum = fabs( pressure.calcSorsum() );
	
	++counter;
    }

    return sorsum;
    
}
